Estimation of vibration amplitude and elastic properties of extra-capillary tissue with ultrasound driven vibration of intra-capillary gas bubbles

ABSTRACT

Estimation of vibration amplitude of intra-capillary micro-bubbles driven to vibrate with an incident ultrasound wave with amplitude and frequency to adjust the drive amplitude of the incident wave to obtain specified vibration amplitude of extra-capillary tissue. Estimation uses transmission of M groups of pulse complexes having low frequency pulse (LF) at bubble drive frequency, and high frequency (HF) pulse with angular frequency ω H &gt;˜5ω L , and pulse duration shorter than π/4ω L  along HF beam. The phase between HF and LF pulses is ω L t m  for each group, where t m  varies between the groups. Within each group, LF pulse varies between pulse complexes in amplitude and/or, where the LF pulse can be zero for a pulse complex, and LF pulse is different from zero for pulse complex within each group. HF receive signals are processed to obtain a parameter relating to bubble vibration amplitude when the HF pulse hits bubble.

RELATED APPLICATION

This application claims priority from U.S. provisional patent application No. 63/248,180 filed on Sep. 24, 2021, the entire content of which is incorporated by reference.

FIELD OF THE INVENTION

The invention addresses methods and instrumentation for using ultrasound vibrations of intra-capillary micro-bubbles to increase fluid to gas phase transition of fluor-carbon (FC) droplets in humans.

DISCUSSION OF RELATED ART

Liquid to gas phase transition of perfluor-carbon nano- and micron-droplets is an emerging field for improving ultrasound diagnosis and assisted therapy of diseases such as cancer, inflammation, atherosclerosis, and myocardial infarction.

Ultrasound stimulated phase transition of ˜3 μm diameter Per-Fluor-Carbon (PFC) droplets produces micro-bubbles of ˜15 μm diameter in larger arteries. In capillaries with typical diameters ˜6-10 μm, such gas bubbles fills the capillary for ˜30-100 μm length with a slow absorption decay into blood over some minutes. Ultrasound vibration of the intra-capillary bubbles produces vibration of the extra-capillary tissue ˜1 μm and has shown to increase efficacy of cancer chemotherapy in pre-clinical studies. There are also indications that such intra-capillary bubble vibrations can increase the immune response to a particular cancer in association with other therapies, and reduce myocardial hypoxia after coronary reopening after a cardiac ischemia attack.

The therapeutic effect depends on adequate amplitude of the capillary wall vibrations. Bubbles obtained by evaporation of PFC droplets of ˜3 μm diameter produce resonance frequencies around ˜500 kHz, and to obtain desirable vibration amplitude of the capillary wall one wants to vibrate the bubbles close to their resonance frequency. The vibration amplitude depends on the amplitude and frequency of the incident ultrasound drive wave, the capillary dimension, and particularly also the shear stiffness of the extra capillary tissue. Especially can both the capillary diameter and the shear stiffness and absorption of the extra-capillary tissue vary by a large amount between different cases. It is therefore desirable to have an online assessment of the vibration amplitude of bubbles that fills the capillary. The invention presents new methods and instrumentation that produces estimates of both the bubble vibration amplitude and its resonance frequency, and also extracts elastic information of tissue surrounding the bubbles.

SUMMARY OF THE INVENTION

An overview of the invention is presented, where the overview is a short form and by no means represents limitations of the invention, which in its broadest aspect is defined by the claims appended hereto.

The invention presents methods and instrumentation for estimation of vibration amplitude of intra-capillary micro-bubbles that fills the capillary for a length and are driven to vibrate with an incident ultrasound wave of given amplitude and frequency. The estimated amplitude is for example used to adjust the drive amplitude of the incident wave to obtain a desired vibration amplitude of the extra-capillary tissue.

The estimation method comprises

-   transmitting towards said micro-bubble a set of limited duration     pulse complexes comprising a low frequency (LF) pulse along a LF     transmit beam (T-beam) with angular frequency at close to the     typical vibration drive frequency, and a high frequency (HF) pulse     with angular frequency ω_(H)>˜5 ω_(L), and pulse duration shorter     than π/4 ω_(L) along a HF transmit beam (T-beam), where said LF     T-beam overlaps the HF T-beam with essentially the same beam     direction, and -   said set of pulse complexes are transmitted in a set of M groups,     each group comprising at least two pulse complexes where the phase     relation between the HF and LF pulses is ω_(L)t_(m) for each group     where t_(m) varies between the groups, and within each group the LF     pulse varies between pulse complexes in one or both of the amplitude     and polarity, where the LF pulse can be zero for a pulse complex,     and the LF pulse is different from zero for at least one pulse     complex, and -   obtaining for each pulse complex within each group, HF receive     signals from the scattered HF waves from said micro-bubble with a HF     receive beam with beam axis either along the HF transmit beam axis     or crossing the HF beam axis, and -   processing for each group said HF receive signals from the     transmitted pulse complexes with HF pulse phase relation ω_(L)t_(m)     the LF pulse, to obtain a parameter relating to the bubble vibration     amplitude when the HF pulse hits the bubble for said each pulse     complex, and -   selecting said parameter from the group that gives the maximal     vibration amplitude as an indication of the vibration amplitude of     said micro-bubble.

The invention further includes instrumentation for carrying through the methods of estimation in a practical situation.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 a Illustrates a spherical micro-bubble with radius a in a larger vessel.

FIG. 1 b Illustrates a micro-bubble filling a part of a capillary used to vibrate the capillary wall and extra capillary tissue

FIG. 2 a Illustrates interface vibration patterns of a capillary surface vibrated by a micro-bubble as shown in FIG. 1 b.

FIG. 2 b Illustrates a numerical simulation of the displacement vibration amplitude as a function of position in the tissue surrounding the vibrating micro-bubble that fill s part of the capillary.

FIG. 3 Illustrates approximation of a 2nd order transfer function Eq. (5, A14) to the transfer function in Eq. (A12) from incident wave amplitude P_(i) at a frequency f to the bubble volume complex vibration amplitude.

FIG. 4 illustrates two typical combined low frequency (LF) and high frequency (HF) ultrasound pulse complexes used to estimate the LF driven vibration amplitude and transfer function parameters of the micro-bubble.

FIG. 5 a-c Illustrates an ultrasound dual frequency array with LF transmit and HF transmit and receive beams typically used for identification of bubble vibration amplitude and tissue elastic properties.

FIG. 6 a Illustrates the variation of H_(b)(K) for a spherical bubble as a function of frequency when the bubble diameter 2a becomes comparable to the incident wave length.

FIG. 6 b-d Illustrates the variation of the H_(b)(K) for a bubble that fills part of a capillary of radius a_(V)=5 μm as a function of angular direction and incident frequency for a cylindrical length 2l_(c).

FIG. 7 a Shows the detection sensitivity function of the average relative radial vibration amplitude for a bubble that fills part of a capillary with relative cylindrical half length l_(c)/a_(V).

FIG. 7 b illustrates the logarithmic amplitude and phase of a 2^(nd) order volume vibration transfer function in Eqs. (5, A14), with illustration of a linear approximation to parts of the phase function.

FIG. 8 a Illustrates the variation of the relative resonance x_(d) of the relative volume transfer function F_(d) in Eq. (A11) as a function of relative cylindrical half length l_(c)/a_(V) and 3 values of the tissue absorption parameter.

FIG. 8 b Illustrates estimated ratio of the tissue shear wave velocity to the capillary radius, c_(s)/a_(V), Eq. (32), as a function of observed bubble resonance frequency f_(d).

FIG. 8 c Illustrates a logarithmic plot of the real component of the tissue shear stiffness as a function of observed bubble resonance frequency f_(d) and 4 capillary radii.

FIG. 9 Illustrates an overview block diagram of a system for carrying through the methods according to the invention.

DETAILED DESCRIPTION OF EMBODIMENTS

This section gives a more detailed description of example embodiments of the invention, where the examples by no means represents limitations of the invention, which in its broadest aspect is defined by the claims appended hereto.

A. Fundamental Bubble Vibration

The invention presents methods and instrumentation for extracting information about ultrasound driven vibrations of micro-bubbles in tissue, and also to use this information for estimation of tissue properties, such as elastic properties, and properties relating to elastic properties.

The invention is relevant for micro-bubbles with diameter less than the capillaries, and also bubbles that in large arteries have a spherical shape with radius a ˜5-30 μm that gives a bubble volume V_(b)=4πα³/3, while in capillaries with radii a_(V)˜3-5 μm the larger bubbles fill at least part of the capillary. Examples of the larger bubbles are shown in FIG. 1 , where 101 shows a spherical bubble in a larger artery 100, while 111 shows the a bubble filling a smaller capillary 110 with radius a_(V) for a cylindrical length 2l_(c) with two semi-spherical endings 112 of radius a_(V).

There is spatially distributed elasticity and mass density of the tissue, the blood, and the gas, that implies wave propagation with reflection at boundaries. The tissue and the blood have both a high volume compression stiffness λ_(t)=1/κ_(t) (κ_(t)—volume compressibility) that upholds pressure waves with propagation velocity c_(ρ)≈1500 m/s for both tissue and blood. The tissue has in addition a low shear/deformation stiffness μ_(t) that produces tissue shear/deformation waves with no volume compression and low propagation velocity c_(s)˜2-20 m/s. The blood and the gas have negligible shear stiffness, and hence no shear waves. The bubble gas has a low volume compression stiffness, but the mass density is also so low that the gas pressure can be approximated as constant across the bubble volume. The bubble gas pressure depends on the bubble volume and its variation as

$\begin{matrix} {P_{g} = {{P_{0}\left( \frac{V_{b}}{V_{b} + {\delta V}} \right)}^{\gamma} \approx {P_{0} - {\gamma P_{0}\frac{\delta V}{V_{b}}}}}} & (1) \end{matrix}$ $\gamma = {\frac{C_{p}}{C_{V}} \approx 1}$ γ = C_(p)/C_(V) ≈ 1 where P₀˜100 kPa is the environmental pressure, V_(b) is the equilibrium bubble volume, δV is the change in bubble volume, C_(p) and C_(V) are the heat capacities at constant pressure and volume. Typical gases have large molecules that makes C_(p)≈C_(V) and γ≈1. The volume-pressure relation is inherently nonlinear for larger changes in the volume, but we shall in our work use the approximate linear variation with δV/V_(b), to study fundamental aspects of the bubble vibration.

The low tissue shear stiffness allows large deformations of the cylindrical bubble surface 111, with large amplitude, short-range shear deformation of the surrounding tissue with comparatively low volume compression. The large vibration amplitude of the close surrounding tissue produces an outward acoustic radiation force (ARF) that has interesting applications for drug transport both across the capillary wall, through the interstitial space between the tissue cells, and across membranes of the cells and other tissue structures. In this context it is important to obtain adequate bubble/tissue vibration amplitude (˜1 μm ) and one target of this invention is to use ultrasound scattering from vibrating bubbles to measure the bubble vibration amplitude produced by an incident ultrasound wave. The invention further uses the measured bubble vibration amplitude to set the amplitude of the incident vibration wave to achieve desired bubble vibration amplitude in the surrounding tissue to achieve the desired therapeutic effect.

A bubble in a large vessel, as 101 in 100, has a practically spherical shape with a spherical vibration pattern. A drive pressure P_(d) that is the difference between the gas pressure P_(g) and the pressure P_(∞) at a point so distant from the bubble that one can neglect the tissue shear vibration, produces a complex volume vibration amplitude □V_(c) of the bubble as [1]

$\begin{matrix} {\frac{\delta V_{c}}{V_{b}} = {{{H_{d}(\omega)}P_{d}{H_{d}(\omega)}} \approx \frac{- K_{d}}{\omega^{2} - \omega_{d}^{2} - {i2\zeta\omega_{d}\omega}}}} & (2) \end{matrix}$ $K_{d} = {{\frac{3}{\rho a^{2}}\omega_{d}^{2}} = {{\frac{4E_{st}}{\rho a^{2}}2{\zeta\omega}_{d}} = \frac{4\eta_{b}}{\rho a^{2}}}}$ where ρ_(b) and η_(b) are the mass density and coefficient of viscosity of the surrounding blood, and E_(st) is the Youngs module of the blood/gas surface layer. We refer to H_(d) as the direct drive transfer function, and ω_(d) as the direct drive angular resonance frequency produced by the interaction between the co-oscillating mass of the surrounding blood and the elasticity of the blood/gas surface layer. Viscous friction in the blood and the surface layer is given by η_(b).

An incident acoustic wave with pressure amplitude P_(i) and angular frequency w produces a variation in the bubble volume and the gas pressure that produces drive pressure amplitude P_(d) as

$\begin{matrix} {P_{d} = {{{{- \gamma}P_{0}\delta V_{c}/V_{b}} - P_{i}} = {{{- \gamma}P_{0}{H_{d}(\omega)}P_{d}} - P_{i}}}} & (3) \end{matrix}$ $P_{d} = \frac{- P_{i}}{1 + {\gamma{H_{d}(\omega)}P_{0}}}$ which gives the complete transfer function H_(V)(ω) from the incident pressure amplitude P_(i) to the complex relative volume amplitude as

$\begin{matrix} {\frac{\delta V_{c}}{V_{b}} = {{- {H_{V}(\omega)}}P_{i}}} & (4) \end{matrix}$ ${H_{V}(\omega)} = \frac{H_{d}(\omega)}{1 + {\gamma{H_{d}(\omega)}P_{0}}}$ With the 2^(nd) order form of the H_(d)(ω) in Eq.(2), we get the following form of H_(V)(ω)

$\begin{matrix} {{H_{V}(\omega)} = \frac{- K_{d}}{\omega^{2} - \omega_{0}^{2} - {i2\zeta\omega_{0}\omega}}} & (5) \end{matrix}$ ω₀² = ω_(d)² + γK_(d)P₀ ζ = ζ_(d)ω_(d)/ω₀

For a bubble that partially fills a capillary with a large cylindrical region as 111 in FIG. 1 b , the vibration pattern is more complex than the spherical vibrations behind the developments in Eqs. (2-5). The difference in material parameters across the vessel wall, i.e. between gas/tissue for the cylindrical bubble surface region (S_(c)), and blood/tissue for the blood surface region (S_(V)), produce interface waves along the vessel wall with propagation velocity c_(i)˜c_(s) in the vascular region and c_(i)˜2c_(s) where the bubble contacts the vessel, as illustrated in FIG. 2 a.

The FIG. 2 a shows simulations of the spatial variation of the vibration amplitude of a bubble that fills part of a capillary (201 diam 2a_(V)=10 μm) together with interface waves along the blood/tissue interface outside the bubble. The upper panel 200 shows the vibration amplitudes at a frequency close to the volume vibration resonance, and the lower panel 210 shows a frequency 6 times the volume resonance frequency. At equilibrium the cylindrical part of the bubble fills the capillary and ends with the semi-spheres 203. Close to the resonance frequency in 200 the bubble surface amplitude 202 is magnified 5 times [5×] while the blood tissue interface wave amplitude 204 is magnified 50 times [50×]. At 6 times the resonance frequency in 210 the bubble surface amplitude 212 is magnified 7.5 times [7.5×] while the blood tissue interface wave 204 is magnified 75 times [75×] interface wave.

FIG. 2 b gives as 221 a grey-scale display of the simulated tissue displacement amplitude around the bubble (220) close to the bubble volume resonance frequency. The volume amplitude can be expressed by the average normal vibration displacement amplitude across the bubble surface S_(b) with area A_(b) as

$\begin{matrix} {{\delta V_{c}} = {A_{b}{\overset{\_}{\psi}}_{n}}} & (6) \end{matrix}$ ${\overset{\_}{\psi}}_{n} = {\frac{1}{A_{b}}{\int\limits_{S_{b}}{d^{2}r_{0}{\psi_{n}\left( {\underline{r}}_{0} \right)}}}}$

Eqs. (12,13,19) shows that for a bubble with dimensions adequately smaller than the incident wavelength the scattered signal is determined by the bubble volume V_(b) and bubble volume compressibility κ_(b). We note from FIG. 2 b that the normal displacement amplitude on the semi-spherical ends is lower than the cylinder radial displacement amplitude, and these two observations stimulates an approximation of the bubble with cylindrical length 2l_(b) and semispherical ends of radius a_(V), with a cylinder of fixed length 2L where the radius a_(V) vibrates. The bubble volume V_(b) for this cylindrical bubble is the same as for the original spherical bubble 101 in a large artery, and the bubble that fills part of a capillary in 111. We have the relation

$\begin{matrix} {V_{b} = {{\frac{4}{3}\pi a^{3}} = {{{\pi a_{v}^{2}2l_{c}} + {\frac{4}{3}\pi a_{v}^{3}}} = {{\pi a_{v}^{2}2LL} = {l_{c} + {\frac{2}{3}a_{v}}}}}}} & (7) \end{matrix}$ $\frac{l_{c}}{a_{v}} = {\frac{2}{3}\left( {\left( \frac{a}{a_{v}} \right)^{3} - 1} \right)}$ $\frac{L}{a_{v}} = {{\frac{l_{c}}{a_{v}} + \frac{2}{3}} = {\frac{2}{3}\left( \frac{a}{a_{v}} \right)^{3}}}$ where 2L is the length of a cylindrical bubble in a capillary with the same volume V_(b). The vibrating area of the cylindrical bubble is

$\begin{matrix} {A_{c} = {{4\pi a_{v}L} = {{{4\pi a_{v}l_{c}} + {\frac{8\pi}{3}a_{v}^{2}}} = {A_{b} - {\frac{4\pi}{3}a_{v}^{2}}}}}} & (8) \end{matrix}$

We notice that Ac is 4πα_(V) ²/3 which is ⅓ of the surface area of the end spheres. However, the effect of this error is reduced because the displacement of the semi-spherical end surfaces is less than the radial vibration of the cylinder surface as shown in FIG. 2 b . We will therefore use the simplifying approximation

$\begin{matrix} {{\delta V_{c}} = {{2\pi a_{v}{\int\limits_{- L}^{L}{{dz} \cdot {\psi_{r}\left( {\underline{r}}_{0} \right)}}}} = {A_{c}{\overset{\_}{\psi}}_{r}}}} & (9) \end{matrix}$ ${\overset{\_}{\psi}}_{r} = {\frac{1}{2L}{\int\limits_{- L}^{L}{{dz} \cdot {\psi_{r}\left( {\underline{r}}_{0} \right)}}}}$ A_(c) = 4πa_(v)L

Analysis of vibration and scattering from the cylindrical bubble is carried through in Appendix A. We note from 210 in FIG. 2 a that at the higher frequency range we get multiple interface wave lengths along the bubble cylindrical surface which complicates the approximation with the 2^(nd) order transfer function in Eq. (5) over a larger frequency range. FIG. 3 a shows as the solid lines a simulation of the magnitude |H_(V) (ω)| (300) and the phase θ_(V) (ω)=∠H_(V) (ω) (301) of H_(V)(□) from Eq. (A10, A12) of Appendix A for l_(c)/a_(v)=5. The dotted lines show the magnitude |Ĥ_(V)(ω)| (302) and the phase {circumflex over (θ)}_(V) (ω) (303) of the 2^(nd) order approximation Ĥ_(V) (ω) of Eq. (A14). The approximation of the 2^(nd) order function to Eq. (A12) is done in the following steps:

-   1. Re{H_(V)(ω₀)}=0 of Eq. (A14) gives the angular resonance     frequency ω₀. -   2. ζ is determined from the differential of the ohase of H_(V) as in     Eq. (25). -   3. ω_(m) and K_(V) is determined from Eq. (27)

In the example in FIG. 3 a we see that f₀≈340 kHz and we can say that the approximation is very acceptable in a range larger than Ω₀=100-1000 kHz=(f₀/3, 3f₀).

To observe the LF vibration amplitude of a bubble driven by the incident LF wave with angular frequency ω_(L), the invention presents methods and instrumentation building on U.S. Pat. Nos. 8,550,998; 9,291,493; 9,939,413; 10,578,725 and 7,727,156; 10,879,867.

An LF pulse shown by example as 401 or 405 in FIG. 4 is typically transmitted along an LF transmit beam (LT-beam) shown by example as 501 of FIGS. 5 a and b , where 502 illustrates a therapy region. We have illustrated the LF transmit beam 501 with a large aperture and focused onto the therapy region 502 which allows higher LF pressure in the therapy region than in front of this region. For adequate therapy we want a certain surface displacement amplitude of the intra-capillary bubbles (FIG. 1 b ). For hard solid tumors, for example with large ingrowth of connective tissue and/or increased intra-tumoral pressure, adequate displacement amplitude requires increased LF pressure, e.g. 1 Mpa which gives an MI ˜1.6 at LF ˜400 kHz. This is seen as a risk MI for an intra-capillary contrast agent bubble with diam ˜3 μm which could get excessive displacement collapse. However for a bubble that fills the capillary, the vibration displacement is more controlled by the surrounding tissue, and for the corresponding spherical bubbles in the larger arteries and arterioles the artery wall is much stronger so that dangerous damage is less likely. In the normal tissue in front of the tumor where the extra-capillary tissue has much lower shear stiffness, one should have lower LF pressure amplitude. The focused LF beam provides such spatial variation of the LF amplitude.

In many situations, for example trans-abdominal and trans-costal applications, the access window to the therapeutic region is limited, and it is desirable to have a HF transmit and receive beam arrangement with the same direction as the LF beam as illustrated in FIG. 5 a , illustrated as 503 for the HF transmit beam and 504 for the HF receive beam. One would then typically use a wide HF transmit beam (503) with well multi-depth focused parallel HF receive beams 504 illustrated by the dots 505. Such multi-depth focusing is today obtained by software receive beam-formers utilizing parallel processing of Graphic Processing Units (GPUs). Hardware receive beamformers with such parallel receive beams also exists to-day.

With adequate ultrasound access windows, it can be useful to use crossing multi/depth focused HF receive beams exemplified by 514 and the dots 515 in FIG. 5 b , with the forward HF transmit beam as 503 in FIG. 5 a . This reduces the length of the scattering wave vector K=k(e_(i) −e _(s)) as defined in Eq. (12). As discussed after Eq.(12) and in relation to FIG. 6 a-d , this length reduction opens for the us higher HF frequencies with controlled reduction in the magnitude of H_(b) and ΔH_(b)/H_(b) value produced by LF variation of the scatterer dimension.

In a particular solution, one transmits a set S_(N) of N≥1 groups G_(n) pulse complexes where one for each group G_(n) transmits at least two pulse complexes comprising a LF manipulation pulse and a HF observation pulse where the LF pulse varies in one of amplitude, and polarity, and phase between each pulse complexes, where the LF pulse might be zero for one pulse complex, and the LF pulse is non-zero for at least one pulse complex. The LF pulses 401 and 405 in FIG. 4 are example LF pulses for two pulse complexes with opposite polarity of the LF pulses for one group G_(n). The position of the HF pulse relative to the LF pulse is for the group G_(n) defined by the timing t_(n) of the HF pulse transmission. The timing varies between each group indicted by the time points . . . , t_(n−1), t_(n), t_(n+1), . . . shown as 404 in FIG. 4 . The envelopes of the LF pulses varies slowly and shown as 404 in FIG. 4 . The HF is typically >˜5*LF. The LF and the HF pulses can have different phase propagation with depth so that the positions t_(n) of the HF pulses shown in FIG. 4 indicates the position of the HF pulses relative to the LF pulse at the depth of the bubble. This change in phase position can for example be obtained by numerical simulations of the combined LF and HF propagation to the bubble, or through measurements in a water or oil tank. For each group G_(n) defined with position t_(n) of the HF pulse, at least two pulse complexes are transmitted with different polarities and/or amplitudes of the LF pulse, where one LF pulse can have zero amplitude and at least one LF pulse has non-zero amplitude. In a preferred embodiment one transmits two pulse complexes with opposite polarities of the LF pulse shown as 401 and 405 in FIG. 4 .

The time between transmit of the HF pulses is determined so that it is larger than the round trip propagation time T_(R) from HF transmit transducer to the micro-bubble and to the HF receive transducer. In particular solutions where one wants long duration LF transmit vibration of the bubble we can use LF pulses with long duration time T_(L), for example longer than NT_(R) where N is the the number of groups G_(n). Several groups G_(n) of HF transmit pulses are then transmitted with time distance larger than T_(R) and t_(n) defining the detailed phase relation between the HF pulse and the LF oscillation of the bubble when the HF pulse hits the bubble.

It is also interesting the use a HF transmit beam 513 that crosses the LF beam 501 in the same direction as the set of HF receive beams 514/515. FIG. 5 c illustrates regions of constant phase relation between the HF and LF pulses in the case of crossing HF and LF transmit beams as in FIG. 5 b , where said regions are along lines normal to the sum of the unit direction vectors of the HF and LF transmit beams. The advantage with the crossing HF transmit beam is that for a fixed phase transmit timing t_(n) between the HF pulse and the phase of the LF oscillation, there will be a spatial variation of said phase relation through the therapy region, and there will be line regions exemplified by 516 in FIG. 5 c along which said phase relation is constant. Said line regions have direction normal to e_(L)+e_(H), where e_(L) is the unit vector along the LF transmit beam axis and e_(H) is the unit vector along the HF transmit beam axis. With adequate density of intra-capillary bubbles in the therapy region 502 there will be intra-capillary bubbles with adequate closeness to these line regions to produce the full amplitude of the bubble volume oscillation and full average bubble surface displacement amplitude ψ as defined in Eq. (9). Selecting the detection signal D_(cm) according to Eqs. (18,20) from these bubbles with maximal D_(cn) provides detection signals with the sought after phase relation between the HF pulses and the bubble vibrations when the HF pulses hits the bubble.

For therapy we generally drive the bubble with the low frequency (LF) pressure wave with angular frequency ω_(L), typically in the neighborhood of the bubble resonance frequency to produce large therapeutic cylindrical radius vibration amplitude (˜1 μm). The LF pulse has an amplitude envelope pP_(Le)(t), where p is a polarity parameter of the LF wave. p=1 implies positive LF amplitude (e.g. 401), p=−1 implies inverted LF amplitude (e.g. 405), while p=0 implies zero LF amplitude. Other values/sequences of the polarity parameter can be used to enhance the signal processing, where the cited patents show examples of p=0, ±1; p=±1, ±2, for example to improve suppression of tissue noise and/or compensate for bubble movements between the HF pulses.

The LF pulse produces a vibration in bubble volume as ζV _(b)(t, ω _(L))=pδV _(Le)(t)cos[ω_(L) t+θ(ω_(L))] δV _(Le)(t)≈|H(ω_(L))|P _(Le)[t+τ _(e)(ω_(L))]≈|H(ω_(L))|P _(Le)(0)  (10)

Where θ(ω_(L))=∠H(ω_(L)) is the phase of H(ω_(L)), and τ_(e) (ω_(L))≈θ′(ω_(L)) is an approximate delay of the envelope produced by the transfer function H(ω) for an adequately narrow band HF pulse. The last approximation of δV_(Le)(t) assumes that he envelope varies slowly relative to the ω_(L) oscillation. To measure the bubble vibration amplitude we transmit HF observation pulses as described in FIG. 4 . The bubble volume vibration when the HF pulse hits the bubble is then pδV _(n) =δV _(b)(t _(n), ω_(L))=pδV _(Le)(t)cos[ω_(L) t _(n)+θ(ω_(L))]  (11)

The scattered pressure field of an incident ultrasound pressure wave with angular frequency ω and amplitude P_(i) from a bubble with volume V_(b) and adequately less extension than the incident wavelength λ, is by the Born approximation obtained as

$\begin{matrix} {{p_{s}\left( {\underline{r},\omega} \right)} \approx {{k^{2}\left( {\frac{\kappa_{b}(\omega)}{\kappa_{t}(\omega)} - {{\underline{e}}_{s}{\underline{e}}_{i}}} \right)}V_{b}{H_{b}\left( \underline{K} \right)}\frac{e^{- {ikr}}}{4\pi r}P_{i}}} & (12) \end{matrix}$ ${\kappa_{b}(\omega)} = {{{- \frac{1}{P_{i}}}\frac{\delta V_{b}}{V_{b}}{\kappa_{t}(\omega)}} = {{- \frac{1}{P_{i}}}\frac{\delta V_{t}}{V_{t}}}}$ ${H_{b}\left( \underline{K} \right)} = {\frac{1}{V_{b}}{\int\limits_{Vb}{d^{3}r_{0}e^{{- i}\underline{K}{\underline{r}}_{0}}}}}$ $\underline{K} = {{{\underline{k}}_{i} - {\underline{k}}_{s}} = {k\left( {{\underline{e}}_{i} - {\underline{e}}_{s}} \right)}}$ $k = {\frac{2\pi}{\lambda} = \frac{\omega}{c_{p}}}$ where

-   r is the radial distance from the bubble center of mass, -   k is the wave vector magnitude of both the incident and scattered     wave, and λ=2π/k is the wave length of both the incident and     scattered waves. -   e _(i) is the unit direction vector of the incident wave     approximated as plane wave across the bubble with amplitude P_(i), -   e _(s) is the unit direction vector of the observed scattered wave, -   c_(p) is the pressure wave propagation velocity, -   δV_(b) is the vibration amplitude of the bubble volume, V_(b), and     δV_(t) is the vibration amplitude of the surrounding tissue volume,     V_(t), both produced by the incident wave, -   κ_(b)(ω) is the bubble volume compressibility, and -   κ_(t)(ω) is the volume compressibility of the surrounding material.

The term κ_(b)/κ_(t) represents scattering due to the deviation of the bubble volume compressibility from the tissue volume compressibility, and the term e _(s) e _(i) represents scattering due the relative difference between the bubble mass density (≈0) and the tissue mass density ρ_(t). In a typical situation |κ_(b)(ω)/κ_(t)|>>|e _(s) e _(i)|, and the term e _(s) e _(i) can be neglected. Using c_(p) ²=1/ρ_(t)κ_(t) we modify Eq. (6) to

$\begin{matrix} {{p_{s}\left( {\underline{r},\omega} \right)} \approx {\rho_{t}\omega^{2}{\kappa_{b}(\omega)}V_{b}{H_{b}\left( \underline{K} \right)}\frac{e^{- {ikr}}}{4\pi r}P_{i}}} & (13) \end{matrix}$ ${\kappa_{b}(\omega)} = {{{- \frac{1}{P_{i}}}\frac{\delta V_{cb}}{V_{b}}} = {H_{V}(\omega)}}$

Both for contrast agent bubbles and bubbles for therapy, we are generally interested in vibrating the bubbles in the low frequency (LF) range Ω₀ around the bubble volume vibration resonance frequency, where the 2^(nd) order transfer functions of Eqs. (5,A14) are useful. To measure the bubble vibration amplitude at ω_(L) we transmit in addition high frequency (HF) observation pulses with angular center frequency ω_(H)˜10 ω_(L)>>ω₀ as illustrated in FIG. 4 . We define this HF range as Ω_(H) The 2^(nd) order transfer function of Eq. (5) is still valid at ω_(H) for the spherical bubble in FIG. 1 a . For the bubbles that fills part of a capillary in FIG. 1 b one gets a multipolar variation of the bubble/tissue interface vibration amplitude ψ_(r)(z, ω), shown as 212 in FIG. 2 a . We therefore use the transfer function H_(V) given in Eqs. (A16,A17) to calculate the bubble volume compressibility κ_(b) at ω_(H)>>ω₀ for these bubbles. This gives different formulas for the spherical and the cylindrical bubbles as

Spherical Bubble:

$\begin{matrix} {{{\kappa_{b}(\omega)} \approx {- \frac{K_{d}}{\omega^{2}}}} = {- \frac{3}{\rho a^{2}\omega^{2}}}} & (14) \end{matrix}$ ${p_{s}\left( {\underline{r},\omega} \right)} \approx {{- {{aH}_{b}\left( \underline{K} \right)}}\frac{e^{- {ikr}}}{r}P_{i}{Y_{s}\left( {a;\omega} \right)}} \sim {{aH}_{b}\left( \underline{K} \right)}$ $\begin{matrix} {{{Cylindrical}{Bubble}:{\kappa}_{b}(\omega)} \approx {- \frac{2{K_{d}\left( {L/a_{v}} \right)}}{\rho a_{v}^{2}\omega^{2}}}} & (15) \end{matrix}$ ${K_{pH}\left( {L/a_{v}} \right)} = {{\frac{L}{a_{v}}{K_{d}\left( {L/a_{v}} \right)}} = {\frac{1}{\pi}{\int{{dq}\frac{{qK}_{1}(q)}{K_{0}(q)}\frac{\sin^{2}qL/a_{v}}{q^{2}}}}}}$ ${p_{s}\left( {\underline{r},\omega} \right)} \approx {{- a_{v}}{K_{pH}\left( {L/a_{v}} \right)}{H_{b}\left( \underline{K} \right)}\frac{e^{- {ikr}}}{r}P_{i}}$ ${Y_{c}\left( {a_{v};\omega} \right)} \sim {a_{v}{K_{pH}\left( {L/a_{v}} \right)}{H_{b}\left( \underline{K} \right)}}$ where it is shown the major dependencies of the receive signal Y (D;ω) on the bubble dimensions D. For the spherical bubble D is the bubble radius a, while for the bubble filling parts of the capillary D is (a_(V),L) or (a_(V), a) by Eq.(7). A change in the bubble dimension/volume at the HF/LF time reference t_(m) when the HF pulse hits the bubble, produces a change in the amplitude of the HF receive signal. For polarity p of the LF pulse the major change in bubble dimensions are Sphere: α=α₀ +pψ Cylindrical: α_(V)=α_(v0) +pψ   (16)

The spatial average ψ of the normal surface displacement is found for the cylinder as described in Eq. (9). The received signal from the HF pulse has angular frequency components in a band Ω_(H) centered around ω_(H). The propagation velocity c_(p) of forward transmitted HF pulse observes a modification by the co-propagating LF wave as described in U.S. Pat. Nos. 9,291,493; 11,280,906; 11,275,006. c _(p) =c ₀(1+pβP _(L))   (17) where c₀ is the low amplitude linear propagation velocity, P_(L) is the pressure of the low frequency wave at the location of the HF pulse, β is a nonlinearity parameter of the tissue bulk compression, and p is a scaling parameter defined in relation to Eqs. (10,11). This nonlinear propagation effect gives a non-linear propagation delay (NPD) to the HF pulse, pτ_(n)(r), and a weaker pulse form distortion (PFD) that increases with the propagation distance r is before scattering, as described in the cited patents. When the transmitted HF pulse is at the crest or through of the LF wave, the NPD has a larger effect than the PFD that generally can be neglected. When the transmitted HF pulse is at a zero crossing of the LF wave, the PFD has a larger effect than the NPD that generally can be neglected.

The back scattered HF signal will generally in addition to the scattering from the bubbles contain considerable scattering from tissue surrounding the bubble, and it is devised to suppress the tissue signal, for example by correcting the received signal for adequate NPD and/or PFD according to the methods in the cited patents, or using the 2^(nd) harmonic component of the received HF signal, or a combination of both. We notice, however that the transmitted LF amplitude to vibrate bubbles is generally lower than for imaging, and both the NPD and the PFD have less effect, and corrections for NPD and/or PFD can in many cases be neglected. Before further processing, we assume that the received HF signals are adequately corrected with the NPD and PFD, as described in the cited patents, and in the notation for the received HF signals below, we assume that that adequate correction have been done.

For illustration of a simplest form, we linearize the receive signals and form a detection signal D_(m) with p=±1 as

$\begin{matrix} {{Y_{np}(\omega)} \approx {{Y_{0}(\omega)} + {p\delta{Y_{n}(\omega)}}}} & (18) \end{matrix}$ ${\delta Y_{n}} = {\frac{\partial Y_{n}}{\partial a_{v}}\psi_{n}}$ ${D_{n}(\omega)} = {\frac{Y_{n -} - Y_{n +}}{Y_{n -} + Y_{n +}} = {\frac{1}{Y_{0}}\frac{\partial Y_{n}}{\partial a_{v}}\psi_{n}}}$ where the subscript n indicates the time point t_(n) in FIG. 4 the measurements are done for, and the subscript n±implies p=±1. Other values for the polarity parameter p can be used as described in relation to FIG. 4 and Eqs. (10,11). We do in most of the equations operate in the angular frequency domain. As understood by any-one skilled in the art, Eq. (18) and other signal processing can through inverse Fourier transform be carried through in the time domain for a large range of depths.

Components for Y_(s) and Y_(c) in Eqs. (14,15) that have low variation with bubble dimension variations are shorted in the ratio of D_(m). In particular has H_(b) low sensitivity to variations in bubble dimensions that are small compared to the wave length. H_(b)(K) is the Fourier transform of the bubble volume with Fourier vector K=k(e _(i)−e _(s)). The largest magnitude K _(max)=2ke _(i)=(4π/λ)e _(i) is found for back-scattering where e _(s)=−e _(i), which the situation in FIG. 5 a where the HF transmit (503) and receive (504) beams have the same. With a HF receive cross beam as 514 in FIG. 5 b and the HF transmit beam as 503 in FIG. 5 a we get K=k|e _(i)−e _(s)|<K_(max)=4π/λ. For a spherical bubble we get

$\begin{matrix}  & (19) \end{matrix}$ $\begin{matrix} {{{Spherical}{bubble}:{H_{b}\left( \underline{K} \right)}} = {\frac{3}{({Ka})^{3}}\left\lbrack {{\sin{Ka}} - {{Ka}\cos{Ka}}} \right\rbrack}} & \left. a \right) \end{matrix}$ $\begin{matrix} {{{Cylindrical}{bubble}:{H_{b}\left( \underline{K} \right)}} = {\frac{2{\sin\left( {{KL}\cos\varphi} \right)}}{{KL}\cos\varphi}\frac{J_{1}\left( {{Ka}_{v}\sin\varphi} \right)}{{Ka}_{v}\sin\varphi}}} & \left. b \right) \end{matrix}$

We notice that for the sphere we have for back scattering Ka=(4πa/c_(p))·f where f is the frequency and c_(p)≈1500 m/s is the propagation velocity of the pressure wave. FIG. 6 a shows ΔH_(b)/H_(b) as a function of frequency for the spherical bubble where the curves 601, . . . , 605 shows the curves for bubble diameters 2a=10, 15, 20, 25, 30 μm, and ψ/a=0.1. ΔH_(b)/H_(b0)<0.02 is accepted as a tolerable error, which gives an upper frequency limit where H_(b) can be accepted as a constant in the processing.

Scattering from the cylindrical bubble in FIG. 1 b gives more complex dependency of H_(b) with bubble dimensions when the bubble long axis is larger than the capillary diameter. H_(b) then becomes a function of both the short axis and long axis bubble dimensions, and the direction of K in relation to the bubble, for example described by the angle φ to the bubble long axis in FIG. 1 b . FIG. 6 b shows H_(b) as a function of φ. The calculations are done for a capillary radius of a_(V)=5 μm, with a bubble cylindrical length 2l_(c)=20 μm in 610, 2l_(c)=30 μm in 610, 2l_(c)=40 μm in 620, 2l_(c)=50 μm in 630. Within each of the circular sub-figures we have a frequency from the outer to the inner curve of f=2.5, 5, 7.5, 10, 15 MHz in each of the Figures. We notice that for K normal to bubble long axis there is negligible change in H_(b) in the actual frequency range, while for K along the bubble long axis the variation is much larger. To approximate H_(b) to constant, we first must both have H_(b) adequately large to get adequate signal in Eq. (14,15,18), while at the same time we must have ΔH_(b)/H_(b0)<0.025 to assume H_(b) constant.

ΔH_(b)/H_(b0) is analyzed in FIG. 6 c and d for the frequencies f=2.5, 5, 7.5, 10, 15 MHz that are displayed from the inner to the outer curve in each circular sub-figure. The capillary radius a_(V)=5 μm in all Figures, and the cylinder half-length l_(c) is 10, 15, 20, 25 μm as listed under each circular sub-figure. FIG. 6 c shows ΔH_(b)/H_(b0) for for Δa_(V)=1 μm, i.e. Δa_(V)/a_(V)=0.2.

For this situation ΔH_(b)/H_(b0)<0.02 is obtained for f=2.5 and 5 MHz. Reducing □a_(V) to 0.5 □m we have ΔH_(b)/H_(b0)<0.02 up to 7.5 MHz for all situations. FIG. 6 d shows ΔH_(b)/H_(b0) for Δl_(c)=1 μm for all values of l_(c). For this situation ΔH_(b)/H_(b0)<0.02 is obtained for f=2.5 and 5 MHz. Reducing Δl_(c) to 0.5 μm we have ΔH_(b)/H_(b0)<0.02 up to 7.5 MHz for all situations.

We expect capillaries and hence also bubbles to have random position within the tumor, and if we select the bubbles with the strongest back-scatter, i.e. incident wave normal to the capilary length and φ close to ±π/2 eg, we can assume that the variation of H_(b) for the variations of the bubble dimension produced by the LF pulse is negligible for all the given bubble lengths and frequencies HF<15 MHz.

We note that when the incident HF beam hits close to the short axis of the bubble, both H_(b0) and ΔH_(β)/H_(b0) are approximately independent of the dimension variations for HF<15 MHz When we can neglect the changes in H_(b) with the bubble dimension variations produced by the polarity change in the incident LF pulse, H_(b) is essentially canceled in the fraction of D_(m) in Eq. (18). Neglecting H_(b) in the expression for Y_(s), Y_(c) in Eq. (14,15) we see that D_(m) of Eq. (18) gives

$\begin{matrix}  & (20) \end{matrix}$ $\begin{matrix} {{{{Spherical}:D_{sm}} = {\frac{Y_{{sm} -} - Y_{{sm} +}}{Y_{{sm} -} + Y_{{sm} +}} = \frac{\psi_{sm}}{a}}},} & \left. a \right) \end{matrix}$ $\begin{matrix} {{{{Cylindrical}:D_{cm}} = {\frac{Y_{m -} - Y_{m +}}{Y_{m -} + Y_{m +}} = {A_{cm}^{- 1}\frac{\psi_{cm}}{a_{v}}}}}{\frac{\psi_{cm}}{a_{v}} = {A_{cm}\frac{Y_{m -} - Y_{m +}}{Y_{m -} + Y_{m +}}}}} & \left. b \right) \end{matrix}$ $A_{cm}^{- 1} \approx {\frac{1}{K_{pH}\left( {L/a_{v}} \right)}\left( {{K_{pH}\left( {L/a_{v}} \right)} + {a_{v}\frac{\partial{K_{pH}\left( {L/a_{v}} \right)}}{\partial a_{v}}}} \right)}$ where K_(pH)(L/α_(V)) is given in Eq. (15, A17). We note that ψ_(cm) is the spatial average of the normal displacement ψ_(r) for the cylinder, as shown in Eq.(A12), while ψ_(sm) is the maximal displacement of the spherical bubble radius, as the bubble displacement is constant across the bubble surface. A_(cm) as a function of l_(c)/a_(V) is shown as 700 in FIG. 7 a . The relation between L and l_(c) is given in Eq. (7). We notice that for l_(c)=0 we have a sphere, and FIG. 7 a shows A_(cm)=1, which is the value given for the sphere in Eq. (21a). This gives trust in the curve in the curve 700 for A_(cm). Newly evaporated ACT bubbles have a cylindrical length l_(c)/a_(V)˜3- 20, depending on capillary radius and degree of gas absorption in the blood. The curve 700 shows that A_(cm)≈0.55-0.65 is a good assessment to obtain estimates of via, with interesting accuracy.

The relative volume variation δV/V_(b) can be approximated from the relative displacement for the spherical and cylindrical bubbles with the linear approximations as

$\begin{matrix} {{{Spherical}:\frac{\delta V_{s}}{V_{b}}} \approx {3\frac{\psi_{sm}}{a}}} & (21) \end{matrix}$ ${{Cylindrical}:\frac{\delta V_{c}}{V_{b}}} \approx {2\frac{\psi_{cm}}{a}}$

Considering measurements errors and relative displacement magnitudes, these approximate expressions have good practical value. More accurate nonlinear expressions between δV/V_(b) and the displacements can be calculated by anyone skilled in the art.

For the 2^(nd) order transfer functions in Eqs. (5,A14) it is possible to measure/estimate the three major parameters, w₀, ζ, and K_(V), for example in the following way. We get a phase estimate θ(ω_(L))=∠H (ω_(L)) of the transfer function at ω_(L) for the t_(n)=t_(m) that gives maximum of ζV_(b) in Eq. (11), i.e. ω_(L) t _(m)+θ(ω_(L))=0⇒θ(ω_(L))=−ω_(L) t _(m)   (22)

FIG. 7 b shows an example amplitude (701) and phase (702) of the 2^(nd) order transfer function H_(V)(ω). Due to the − sign in front of H_(V)(ω) in Eqs. (5,A14) the resonance frequency ω₀ the phase of the volume vibration amplitude δV_(cb)(ω)is −π/2 behind the driving pressure amplitude −P_(i), while for ω adequately less than ω₀ the phase-lag is ≈0 behind −P_(i), i.e. the vibration amplitude has close to opposite polarity of P_(i). For ω adequately larger than ω₀, the phase-lag is ≈−π behind −P_(i), i.e. the δV_(cb) (ω) has the same polarity as P_(i). The amplitude and phase of the bubble vibration hence depend on the details of the transfer function (Gain K, resonance frequency ω₀ and absorption ζ).

Measuring {circumflex over (θ)}(ω_(L1)) and {circumflex over (θ)}(ω_(L2)) from Eq. (22) for two low frequencies ω_(L1) and ω_(L2) where we get the maximal detection signal for t_(m1) and t_(m2), we can for example though linear interpolation, indicated as 703 in FIG. 7 b , write the phase variation between ω_(L1) and ω_(L2) as

$\begin{matrix} {{\hat{\theta}(\omega)} = {{\hat{\theta}\left( \omega_{L1} \right)} + {\frac{{\hat{\theta}\left( \omega_{L2} \right)} - {\hat{\theta}\left( \omega_{L1} \right)}}{\omega_{L2} - \omega_{L1}}\left( {\omega - \omega_{L1}} \right)}}} & (23) \end{matrix}$

At the resonance frequency □₀ we shall have θ(ω₀)=−π/2 which gives an estimated resonance frequency from Eq. (23) as

$\begin{matrix} {{\hat{\omega}}_{0} = {\omega_{L1} - {\frac{{\pi/2} + {\hat{\theta}\left( \omega_{L1} \right)}}{{\hat{\theta}\left( \omega_{L2} \right)} - {\hat{\theta}\left( \omega_{L1} \right)}}\left( {\omega_{L2} - \omega_{L1}} \right)}}} & (24) \end{matrix}$ Matching the differential of the phase in Eq. (5,A14) at ω=ω₀ to Eq. (24) we get an estimate for the absorption losses in the tissue vibrations as

$\begin{matrix} {\zeta = {\left. {- \frac{1}{\omega_{0}{\theta^{\prime}\left( \omega_{0} \right)}}}\Rightarrow\hat{\zeta} \right. = \frac{\omega_{L2} - \omega_{L1}}{{\hat{\omega}}_{0}\left( {{\theta\left( \omega_{L2} \right)} - {\theta\left( \omega_{L1} \right)}} \right)}}} & (25) \end{matrix}$ Another estimate of the absorption factor ζ is from Eq. (5) also found for a single low frequency ω_(L) as

$\begin{matrix} {\hat{\zeta} = {\frac{1}{2}\left( {\frac{\omega_{L}}{{\hat{\omega}}_{0}} - \frac{{\hat{\omega}}_{0}}{\omega_{L}}} \right)\tan{\hat{\theta}\left( \omega_{L} \right)}}} & (26) \end{matrix}$ The parameter K_(V) can for example be found by the angular frequency ω_(m) that gives maximal amplitude of H_(V)(ω). From the estimates of {circumflex over (ω)}₀ and {circumflex over (ζ)} and we can obtain an estimate of {circumflex over (ω)}_(m) ²={circumflex over (ω)}₀ ²(1−2{circumflex over (ζ)}²) {circumflex over (K)} _(V) =|H _(V)({circumflex over (ω)}_(m))|√{square root over ({circumflex over (ω)}₀ ⁴−{circumflex over (ω)}_(m) ⁴)}  (27) where we measure |H_(V)| at the estimated angular frequency ω_(m). If we do not have an estimate of {circumflex over (ζ)}, we can search for the c that gives max |H_(V)| and do the following calculations

$\begin{matrix} {{\hat{K}}_{V} = {{{❘{H_{V}\left( {\hat{\omega}}_{m} \right)}❘}\sqrt{{\hat{\omega}}_{0}^{4} - {\hat{\omega}}_{m}^{4}}\hat{\zeta}} = \sqrt{\frac{\omega_{0}^{2} - \omega_{m}^{2}}{2\omega_{0}^{2}}}}} & (28) \end{matrix}$

Estimating the phases {circumflex over (θ)}_(n)={circumflex over (θ)}(ω_(Ln)) at more than two frequencies ω_(Ln), n=1, 2, . . . , N, we can improve the estimate of {circumflex over (ω)}₀ with more complex interpolations, as known by any-one skilled in the art. By example, for n=1, . . . , N we the set of measured frequencies by the vector ω _(L)={ω_(L1), . . . , ω_(Ln), . . . , ω_(LN)} and the set of estimated phases by ω _(L)={ω_(L1), . . . , ω_(Ln), . . . , ω_(LN)}, we can then form an estimate of {circumflex over (ω)}₀ and {circumflex over (ζ)}{circumflex over (ω)}₀ through minimization of a mean square functional as

$\begin{matrix} {{J\left( {\omega_{0},{\zeta;{\underline{\omega}}_{L}}} \right)} = {\left\lbrack {{\underline{\hat{\theta}}\left( {\underline{\omega}}_{L} \right)} - {\underline{\theta}\left( {\omega_{0};\zeta;{\underline{\omega}}_{L}} \right)}} \right\rbrack^{T}\left\lbrack {{\underline{\hat{\theta}}\left( {\underline{\omega}}_{L} \right)} - {\underline{\theta}\left( {\omega_{0},{\zeta;{\underline{\omega}}_{L}}} \right)}} \right\rbrack}} & (29) \end{matrix}$ $\left( {{\hat{\omega}}_{0},\hat{\zeta}} \right) = {\min\limits_{\omega_{0},\zeta}{J\left( {\omega_{0},{\zeta;{\underline{\omega}}_{L}}} \right)}}$

Other forms of functionals for minimization, also including added information, can be formed by any-one skilled in the art.

From Eq. (5-7) we notice that for {circumflex over (θ)}(ω_(L))+ω_(L) _(t) _(m)≈0 the amplitude of the received HF signal takes the form Y _(np) =Y _(n0) +pδY _(n) ˜p|H _(V)(ω_(Ln))|P _(Le)(t _(n))   (30) and we can for example include the amplitude of the received signal ˜|H_(V)(ω_(Ln))| for the different {ω_(L1), . . . , ω_(LN)} in addition to the phase in the functional for parameter estimation in Eq. (29).

The extra capillary tissue is determined by the mass density, ρ_(t), and the real and imaginary part, μ_(tr) and μ_(ti) of the complex shear stiffness of the tissue, as shown in Eq. (A6, A7). We can approximate ρ_(t)≈1100 kg/m³ for most solid tumor tissues, and hence it is μ_(tr) and μ_(ti) which are the most sought after tissue parameters. From the measured/estimated {circumflex over (ω)}₀ and {circumflex over (K)}_(V), we see from Eq. (A18) that we can obtain an estimate of {circumflex over (ω)}_(d) as

$\begin{matrix} \begin{matrix} {{\hat{\omega}}_{d} = \sqrt{{\hat{\omega}}_{0}^{2} - {\gamma{\hat{K}}_{V}P_{0}}}} & {\zeta_{d} = {\frac{\omega_{0}}{\omega_{d}}\zeta}} & {{\hat{x}}_{d} = \frac{{\hat{\omega}}_{d}}{\omega_{s}}} & {\omega_{s} = {\frac{c_{s}}{a_{v}} = {\frac{1}{a_{v}}\sqrt{\frac{\mu_{tr}}{\rho}}}}} \end{matrix} & (31) \end{matrix}$ where γ≈1 and P₀≈101 kPa.

From F_(d)(x) in Eq. (A11) we obtain estimated x_(d)(l_(c)/a_(V)) from Re{F_(d)(x_(d))=0. FIG. 8 a shows x_(d)(l_(c)/a_(V)) for a set of absorption parameters α=0.5, 1.0, 1.5, defined from η=μ_(ti)/μ_(tr)=2αx_(d). We notice that the variation of the curves with α Is barely visible, only for l_(c)/a_(V)<4. We also notice from Eq. (A10) that ω_(s)=c_(s)/a_(V)=ω_(d)/x_(d) and μ_(tr)=ρα_(V) ²ω_(s) ² which allows us to calculate

$\begin{matrix} \begin{matrix} {{\frac{c_{s}}{a_{v}}\left( {l_{c}/a_{v}} \right)} = {\omega_{s} = \frac{\omega_{d}}{{\hat{x}}_{d}\left( {l_{c}/a_{v}} \right)}}} & {{\mu_{tr}\left( \omega_{d} \right)} = {\rho a_{v}^{2}\frac{\omega_{d}^{2}}{{\hat{x}}_{d}^{2}\left( {l_{c}/a_{v}} \right)}}} \end{matrix} & (32) \end{matrix}$

FIG. 8 b shows ω=c_(s)/a_(V) as a function of {circumflex over (ω)}_(d) for of l_(c)/a_(V)=0, 4, and 10 and FIG. 8 c shows μ_(tr)(ω_(d)) for a_(V)=2, 4, 5, and 6 μm.

The power dissipated by the surface vibration at a surface area element d²r is

$\begin{matrix} {{dW} = {{{- \frac{1}{2}}{Re}\left\{ {d^{2}{ri}\omega{\psi_{n}\left( \underline{r} \right)}^{*}P_{i}} \right\}} = {{- \frac{\omega}{2}}{Im}\left\{ {d^{2}r{\psi_{n}\left( \underline{r} \right)}P_{i}} \right\}}}} & (33) \end{matrix}$ We assume P_(i) approximately constant across the bubble surface, and integrated across the whole bubble surface gives the following power dissipation that provides a heating energy per unit time

$\begin{matrix} {W = {{{- \frac{\omega}{2}}{Im}\left\{ {P_{i}{\int{d^{2}r{\psi_{n}\left( \underline{r} \right)}}}} \right\}} = {{- \frac{\omega}{2}}{Im}\left\{ {P_{i}\delta V_{cb}} \right\}}}} & (34) \end{matrix}$ Inserting the transfer function from Eq. (5), we get

$\begin{matrix} {{W(\omega)} = {{\frac{\omega}{2}V_{b}P_{i}^{2}{Im}\left\{ {H_{V}(\omega)} \right\}} = {\frac{\omega_{0}\zeta K_{V}V_{b}\omega^{2}/\omega_{0}^{2}}{\left( {1 - {\omega^{2}/\omega_{0}^{2}}} \right) + {4\zeta^{2}\omega^{2}/\omega_{0}^{2}}}P_{i}^{2}}}} & (35) \end{matrix}$

Hence, from estimation of ω₀ and ζ and assessment of the bubble volume and the incident pressure amplitude, we can assess the dissipated power from the bubble on the surrounding tissue at ω=ω_(L). Because ω_(L) is low, we can estimate the low frequency amplitude P_(i)=P_(L) using an assessment of the low absorption. From the ultrasound image (2D or 3D) we can automatically estimate the number of bubbles N per unit volume and assess the total heat generation as W _(T)(ω_(L))=NW(ω_(L))   (36) that can be used to estimate temperature increase in the tissue to avoid uncontrollable results.

FIG. 9 shows a block diagram of an example instrument for carrying out the methods according to the invention. 900 shows by example a 3D ultrasound probe comprising a dual frequency linear array 901 with a set of M LF elements and N HF elements in an azimuth direction indicated by the arrows 906. The dual frequency band linear array can be made according to known methods, for example as described in U.S. Pat. Nos. 7,727,156; 10,879,867. The LF and HF elements of the array are via a cable 902 connected to a transmit/receive unit 903 that connects each LF array element to LF transmit circuits, and each HF element to HF transmit/receive circuits comprising a HF transmit amplifier and a HF receive amplifier where the output of the HF receive amplifier is further connected to an analog to digital converter (A/D) presenting a digital representation of the HF received signals from all HF receive elements, according to known methods. The AD converter can in a modified embodiment present digital representations of the I-Q components of the HF receive signals from each HF element that represents the same information as the radio frequency (RF) HF signal, according to known methods. In a modified version of the instrument where one wants to receive the scattered LF pulse from the bubble, the LF elements are also connected to LF transmit/receive circuits.

For 3D scanning of the ultrasound beams, the linear array 901 can in this example embodiment be rotated around the long axis 904 that provides a mechanical scanning of the LF/HF beam in an elevation direction, indicated by the arrows 905. For each elevation position of the array, one does electronic scanning of LF/HF transmit beams in an azimuth direction indicated by the arrows 906, through electronic selection of transmitting LF and HF elements, and transmitting combined LF/HF pulse complexes for example as shown in FIG. 4 , with selected beam directions and foci. An example HF transmit beam is illustrated schematically as 907 within a 2D azimuth plane 908 with given elevation position within a total 3D scan volume 909. Alternative elevation movements of the array, like side-ways movement can equivalently be done according to known methods, depending on the space available for such movement, and the shape of the object. Pure 2D scanning of the beam, and even measurements along a single beam direction is also within the scope of the invention. The array might also have a 2D matrix distribution of elements for 3D focusing and scanning of the beams, as known to anyone skilled in the art.

At least two pulse complexes with different LF pulses, for example as illustrated in FIG. 4 , are transmitted for each transmit beam direction. The LF pulse might be zero in one pulse complex per HF transmit beam, but must be non/zero in at least one pulse complex for each HF transmit beam direction.

Two versions of the instrument are useful, where the first version 903 comprises beam former for HF receive cross-beams, illustrated as 914 in the 2D scan plane 908, and HF back scatter receive beams with the same axis as the HF transmit beam 907, shown in further detail in FIG. 5 . In a preferred embodiment the HF back-scatter receive beam is equal to the HF transmit beam as this improves suppression of multiple scattering noise in the HF back-scatter receive signal, as discussed in U.S. Pat. Nos. 8,550,998; 9,291,493; 9,939,413; 10,578,725, and synthetic depth optimized focus is obtained through lateral filtering of the HF receive signals. During the scan, the HF cross-beam and back-scatter receive signals are via the high speed bus 910 transferred to the processor 911 for storage and further processing.

The processor 911 typically comprises at least one multicore central processing unit (CPU) and at least one graphics processor unit (GPU) that are SW programmable, according to known methods. The processor receives user inputs from a user/operator input unit 913 that operates according to known methods, and displays image data and other information necessary for communication with the user/operator through a combined display and audio unit 912, according to known methods.

In the second version of the instrument, the digital HF receive signals from each HF receive elements and each transmitted pulse complex are via the high-speed bus 910 transferred to the processor 911 for storage and further processing. For 2D imaging in the second version, a SW program in the processor unit 911 combines HF receive signals from multiple HF receive elements to produce both HF backscatter receive beams and HF receive cross-beams crossing each HF transmit beam in the 2D set, for example as shown as 914 and in more detail in FIGS. 5 b and c . A set of HF back-scatter receive signals from HF back-scatter receive beams with the same axis as the HF transmit beams, and preferably also equal to the HF transmit beams 907, also shown in FIG. 5 , can be produced by a software program.

One use of the cross-beam (914) HF receive signals is to estimate the nonlinear propagation delay at a multitude of depths along the HF transmit beam with low influence of multiple scattered noise, as described in US Pat Appl 2020/0405268. A set of cross beams is also useful to provide better suppression of signal from tissue, that is noise for detecting bubbles, vibration amplitudes, and tissue parameters. When we have imaging situations where the level of multiple scattering noise in the back-scatter HF receive signal is low compared to the 1st order scattered signals, the nonlinear propagation delay at multitude of depths along the HF transmit beam (907) can be obtained from the back-scatter HF receive signal, as described in the above cited US Patents, and the beam former for HF receive cross-beams, has less importance and can be omitted.

The processor 911 further comprises means, for example implemented as SW programs, to carry through the processing as described above to estimate

-   i) the relative variation of the bubble volume or displacement     amplitude when driven by the LF signal, as for example described in     Eqs. (19,21) and following equations, and -   ii) the bubble vibration resonance frequency and absorption, as for     example described in Eqs. (22-25, 30) and related equations, and -   iii) parameters of shear elasticity of the extravascular tissue, for     example as shown in FIGS. 8 b and c.

The method and instrument is useful for estimating vibration amplitudes of bubbles of any dimension at an angular frequency we define as the low frequency, ω_(L). For spherical bubbles we have a good analytic relationship of a 2^(nd) order transfer function, Eq. (5) between the bubble radius and the scattered wave over a large frequency range from well below to well above the resonance frequency. For Bubbles that fill a capillary, like in FIG. 1 b , we have in Appendix A develop an analytic expression for H_(V)(w) in Eq. (A12) with an approximate 2^(nd) order transfer function for frequencies close to the volume resonance frequency in Eq. (A14). Eqs. (22-32) gives methods to identify the parameters of said model and the tissue surrounding a bubble that fills part of a capillary.

Thus, while there have shown and described and pointed out fundamental novel features of the invention as applied to preferred embodiments thereof, it will be understood that various omissions and substitutions and changes in the form and details of the devices illustrated, and in their operation, may be made by those skilled in the art without departing from the spirit of the invention.

It is also expressly intended that all combinations of those elements and/or method steps which perform substantially the same function in substantially the same way to achieve the same results, are within the scope of the invention. Moreover, it should be recognized that structures and/or elements and/or method steps shown and/or described in connection with any disclosed form or embodiment of the invention may be incorporated in any other disclosed or described or suggested form or embodiment as a general matter of design choice. It is the intention, therefore, to be limited only as indicated by the scope of the claims appended hereto.

APPENDIX A

For the vibration field outside the cylindrical capillary from the vibrating bubble in FIG. 1 b , we use the cylindrical (r,z) coordinates as shown in the Figure. The incompressible field in the vabrating co-oscillating mass outside the capillary can for vibration of a cylindrical bubble with length 2L at an angular frequency ω, be approximated by the Laplace equation for a potential A(r,z;ω)

$\begin{matrix} {{{{\nabla^{2}{A\left( {r,{z;\omega}} \right)}} + {k^{2}{A\left( {r,{z;\omega}} \right)}}}\underset{\lambda\operatorname{>>}{2L}}{\rightarrow}{\nabla^{2}{A\left( {r,{z;\omega}} \right)}}} = 0} & ({A1}) \end{matrix}$ where λ=2πc_(p)/ω is the wavelength of the pressure wave with propagation velocity c_(p). The particle vibration displacement ψ and the pressure P in the co-oscillating mass are obtained from A as ψ(r,z;ω)=∇A(r,z;ω) P(r,z;ω)=ρω² A(r,zω)   (A2) Fourier transform along the z-axis gives

$\begin{matrix} {{{{\frac{1}{r}{\frac{\partial}{\partial r}\left( {r\frac{\partial{A\left( {r,z} \right)}}{\partial r}} \right)}} + \frac{\partial^{2}{A\left( {r,z} \right)}}{\partial z^{2}}}\underset{{FT}_{z}}{\rightarrow}{\frac{\partial^{2}{A\left( {r,k_{z}} \right)}}{\partial r^{2}} + {\frac{1}{r}\frac{\partial{A\left( {r,k_{z}} \right)}}{\partial r}} - {k_{z}^{2}{A\left( {r,k_{z}} \right)}}}} = 0} & ({A3}) \end{matrix}$ Setting x=k_(z)r and A(r,k_(z))=y(k_(z)r) we get the modified Bessel equation with the solution K₀(x) x ² y″(x)+xy′(x)−x ² y(x)=0 x=k _(z) r y(x)=K ₀(x)   (A4) With A defined for r=a_(V) and a tissue mass density ρ, we get the following field amplitudes in the tissue outside the capillary

$\begin{matrix} {{A\left( {r,k_{z}} \right)} = {\frac{A_{0}\left( {a_{v},k_{z}} \right)}{K_{0}\left( {k_{z}a_{v}} \right)}{K_{0}\left( {k_{z}r} \right)}}} & ({A5}) \end{matrix}$ P(r, k_(z)) = ρω²A(r, k_(z)) ${\psi_{r}\left( {r,k_{z}} \right)} = {\frac{\partial{A\left( {r,k_{z}} \right)}}{\partial r} = {{- k_{z}}\frac{A\left( {a_{v},k_{z}} \right)}{K_{0}\left( {k_{z}a_{v}} \right)}{K_{1}\left( {k_{z}r} \right)}}}$ ${\psi_{z}\left( {r,k_{z}} \right)} = {{{- {ik}_{z}}{A\left( {r,k_{z}} \right)}} = {{- {ik}_{z}}\frac{A\left( {a_{v},k_{z}} \right)}{K_{0}\left( {k_{z}a_{v}} \right)}{K_{0}\left( {k_{z}r} \right)}}}$ Inside the co-oscillating mass the incompressible tissue deformation produces shear stresses, which balances to zero inside the tissue, but at the boarder to the bubble at the capillary wall we get a stress contribution to the tissue surface as

$\begin{matrix} {{\sigma_{\mu r}\left( {a_{v},k_{z}} \right)} = {{2\mu_{t}\frac{\partial\psi_{r}}{\partial r}} = {{2\mu_{t}\frac{\partial^{2}A}{\partial r^{2}}} = {{- 2}{\mu_{t}\left( {{\frac{1}{r}\frac{\partial A}{\partial r}} - {k_{z}^{2}A}} \right)}}}}} & ({A6}) \end{matrix}$ ${\sigma_{\mu r}\left( {a_{v},k_{z}} \right)} = {2{\mu_{t}\left( {k_{z}^{2} + \frac{k_{z}{K_{1}\left( {k_{z}a_{v}} \right)}}{a_{v}{K_{0}\left( {k_{z}a_{v}} \right)}}} \right)}{A\left( {a_{v},k_{z}} \right)}}$ where μ_(t)=μ_(tr)+iμ_(ti)=μ_(tr)(1+iη) is the complex shear stiffness of the tissue where the imaginary component produces absorption of oscillation power. At the bubble-tissue surface interface this stress adds to the pressure in Eq. (A5) to the following pressure from the cylindrical bubble surface towards the tissue

$\begin{matrix} {{P_{d}\left( {a_{v},k_{z}} \right)} = {{{\rho{\omega^{2}\left( {a_{v},k_{z}} \right)}} - {\sigma_{\mu r}\left( {a_{v},k_{z}} \right)}} = {\rho{H_{\omega}\left( {\omega;{k_{z}a_{v}}} \right)}{A\left( {a_{v},k_{z}} \right)}}}} & {A(7)} \end{matrix}$ ${H_{\omega}\left( {{\omega;k_{z}},a_{v}} \right)} = {\omega^{2} - {2{\omega_{c}^{2}\left( {\left( {k_{z}a_{v}} \right)^{2} + \frac{k_{z}a_{v}{K_{1}\left( {k_{z}a_{v}} \right)}}{K_{0}\left( {k_{z}a_{v}} \right)}} \right)}}}$ $\omega_{c}^{2} = {{\frac{\mu_{z}}{\rho a_{v}^{2}}\mu_{t}} = {{{\mu_{tr}\left( {1 - {i\eta}} \right)}\omega_{c}^{2}} = {{{\omega_{s}^{2}\left( {1 = {i\eta}} \right)}\omega_{s}^{2}} = \frac{\mu_{tr}}{\rho a_{v}^{2}}}}}$

When the bubble cylindrical length 2l_(c) is adequately larger than the capillary radius a_(V), we can approximate the bubble with semi-spherical end surfaces with radius a_(V) with a cylinder with length 2L as shown in Eqs. (7, 8) above. P_(d)(a_(V), k_(z)) is the z-Fourier transform of a spatially constant drive pressure amplitude P_(d) relative to infinity in the interval (−L, L), and zero outside

$\begin{matrix} {{P_{d}\left( {a_{v},k_{z}} \right)} = {2\frac{\sin k_{z}L}{k_{z}}P_{d}}} & ({A8}) \end{matrix}$ ${A\left( {a_{v},k_{z}} \right)} = {\frac{P_{d}\left( {a_{v},k_{z}} \right)}{\rho{H_{\omega}\left( {\omega;{k_{z}a_{v}}} \right)}} = {\frac{2}{\rho{H_{\omega}\left( {\omega;{k_{z}a_{v}}} \right)}}\frac{\sin k_{z}L}{k_{z}}P_{d}}}$ ${\psi_{r}\left( {a_{v},k_{z}} \right)} = {{{- \frac{k_{z}{K_{1}\left( {k_{z}a_{v}} \right)}}{K_{0}\left( {k_{z}a_{v}} \right)}}{A\left( {a_{v},k_{z}} \right)}} = {{- \frac{k_{z}{K_{1}\left( {k_{z}a_{v}} \right)}}{K_{0}\left( {k_{z}a_{v}} \right)}}\frac{2}{\rho{H_{\omega}\left( {\omega;{k_{z}a_{v}}} \right)}}\frac{\sin k_{z}L}{k_{z}}P_{d}}}$ ${\psi_{r}\left( {a_{v},k_{z}} \right)} = {{- \frac{1}{2\pi}}{\int{{dk}_{z}\frac{k_{z}{K_{1}\left( {k_{z}a_{v}} \right)}}{K_{0}\left( {k_{z}a_{v}} \right)}\frac{2}{\rho{H_{\omega}\left( {\omega;{k_{z}a_{v}}} \right)}}\frac{\sin k_{z}L}{k_{z}}e^{{- {ik}_{z}}z}P_{d}}}}$ The ψ vibration amplitude at the cylindrical part of the bubble amplitude, gives a volume change of the cylinder

$\begin{matrix} {\frac{\delta V_{c}}{V_{b}} = {{\frac{4\pi a_{v}L}{2\pi a_{v}^{2}L}\frac{1}{2L}{\int\limits_{- L}^{L}{{dz}{\psi_{r}\left( {a_{v},z} \right)}}}} = {{- \frac{2}{\rho a_{v}^{2}}}\frac{1}{\pi}{\int{{dk}_{z}\frac{k_{z}{K_{1}\left( {k_{z}a_{v}} \right)}}{K_{0}\left( {k_{z}a_{v}} \right)}\frac{1}{H_{\omega}\left( {\omega;{k_{z}a_{v}}} \right)}\frac{\sin k_{z}L}{k_{z}^{2}L}P_{d}}}}}} & ({A9}) \end{matrix}$ Changing integration coordinates q=k_(z)α_(V), k_(z)=q/α_(V), dk_(z)=dq/α_(V) we get the direct drive transfer function H_(d) as

$\begin{matrix} {\frac{\delta V_{c}}{V_{b}} = {{H_{d}\left( {\omega;{L/a_{v}}} \right)}\frac{P_{d}}{\mu_{tr}}}} & ({A10}) \end{matrix}$ ${H_{d}\left( {\omega;{L/a_{v}}} \right)} = {{- \frac{1}{\pi}}{\int{{dq}\frac{{qK}_{1}(q)}{K_{0}(q)}\frac{2\omega_{s}^{2}}{H_{\omega}\left( {\omega;q} \right)}\frac{\sin^{2}{{qL}/a_{v}}}{q^{2}{L/a_{v}}}}}}$ $\omega_{s}^{2} = {\frac{\mu_{tr}}{\rho a_{v}^{2}} = {{\frac{c_{s}^{2}}{a_{v}^{2}}c_{s}^{2}} = \frac{\mu_{tr}}{\rho}}}$ where c_(s) is the shear wave velocity and ω_(s) is the shear wave angular frequency. Introducing x=ω/ω_(s) we obtain the scaled direct drive transfer function F_(d) as

$\begin{matrix} {\frac{\delta V_{c}}{V_{b}} = {{F_{d}\left( {x;{L/a_{v}}} \right)}\frac{P_{d}}{\mu_{tr}}}} & ({A11}) \end{matrix}$ ${F_{d}\left( {x;{L/a_{v}}} \right)} = {{- \frac{1}{\pi}}{\int{{dq}\frac{{qK}_{1}(q)}{K_{0}(q)}\frac{2}{H_{x}\left( {x;q} \right)}\frac{\sin^{2}{{qL}/a_{v}}}{q^{2}{L/a_{v}}}}}}$ H_(x)(x; q) = x² − 2[q² + qK₁(q)/K₀(q)](1 + iη)

Note that for the cylindrical bubble in a capillary we have above put the surrounding tissue real shear stiffness μ_(tr) outside H_(d) and F_(d). This gives F_(d) the nice property that it does not depend on μ_(tr) and ω_(s), only on η. We assume a linear variation η with frequency, i.e. η=2αx. The resonance frequency of F_(d) is defined as the value x_(d) that equal zero, i.e. Re{F_(d)(x_(d);L/a_(V)}. FIG. 3 b shows x_(d)(l_(c)/a_(V)) as a gives the real part of F_(d) function of l_(c)/a_(V). The angular resonance frequency of H_(d) is ω_(d)=ω_(s)x_(d). The variation of x_(d) with l_(c)/a_(V) is quite low and we can approximate x_(d)≈1 for l_(c)/a_(V) in the range 1 to 8. This allows us to approximate the direct drive resonance function ω_(d)≈ω_(s). For l_(c)/a_(V)=0 we get the direct drive values for the spherical bubble that touches the capillary wall at one circle as x_(d)≈1.7-2 for the absorption α=0.5-1.5 and ω_(d)≈(1.7-2)ω_(s). The resonance frequency for a spherical bubble that is fully confined in soft tissue is ω_(d)≈2ω_(s). The value for a spherical confined in blood with the short touch to the capillary wall does not seem a bad estimate.

Including gas elasticity according to Eq. (3) we get the complete transfer function from incident pressure to relative volume amplitude

$\begin{matrix} {\frac{\delta V_{c}}{V_{b}} = {{{H_{d}\left( {\omega;{L/a_{v}}} \right)}\frac{P_{d}}{\mu_{tr}}} = {{- {H_{v}\left( {\omega;{L/a_{v}}} \right)}}P_{i}}}} & ({A12}) \end{matrix}$ ${H_{v}\left( {\omega;{L/a_{v}}} \right)} = \frac{{H_{d}\left( {\omega;{L/a_{v}}} \right)}/\mu_{tr}}{1 - {\gamma{H_{d}\left( {\omega;{L/a_{v}}} \right)}{P_{0}/\mu_{tr}}}}$

FIG. 3 shows that we for a practical interval around the resonance frequency can approximate the transfer function with a 2^(nd) order transfer function with ω_(d)=ω_(s)X_(d)

$\begin{matrix} {{{\hat{F}}_{d}\left( {x;{L/a_{v}}} \right)} \approx \frac{- K_{F}}{x^{2} - x_{d}^{2} - {i2\zeta_{d}x_{d}x}}} & ({A13}) \end{matrix}$ ${{\hat{H}}_{d}\left( {\omega;{L/a_{v}}} \right)} \approx \frac{- K_{d}}{\omega^{2} - \omega_{d}^{2} - {i2\zeta_{d}\omega_{d}\omega}}$ ω_(d)² = ω_(s)²x_(d)²K_(d) = ω_(s)²K_(F)

Due to the frequency variation of the displacement along the bubble surface as shown in 212 of FIG. 2 a , K_(F) has to be adapted numerically, for example so that the maximal amplitude of the F_(d) approximation in Eq. (A13) matches to the maximal amplitude in Eq. (A11). From Eq. (A12) we see that with a given incident amplitude P_(i) we can include the variations in the gas pressure on the bubble volume to get the complete 2^(nd) order transfer function as

$\begin{matrix} {{{\hat{H}}_{V}\left( {\omega;{L/a_{v}}} \right)} = {\frac{{{\hat{H}}_{d}\left( {\omega;{L/a_{v}}} \right)}/\mu_{tr}}{1 - {\gamma{{\hat{H}}_{d}\left( {\omega;{L/a_{v}}} \right)}{P_{0}/\mu_{tr}}}} = \frac{- K_{V}}{\omega^{2} - \omega_{0}^{2} - {i2{\zeta\omega}_{0}\omega}}}} & ({A14}) \end{matrix}$ $K_{V} = {\frac{K_{d}}{\mu_{tr}} = \frac{K_{F}}{\rho a_{v}^{2}}}$ ω₀² = ω_(d)² + γK_(v)P₀ = ω_(s)²(x_(d)² + K_(F)(L/a_(v))P₀/μ_(tr)) ζ = ζ_(d)ω/ω₀K_(V) = ω_(s)²K_(F)/μ_(tr)

With the 2^(nd) order approximation of H_(V) in Eq. (A14) we can from Eq. (13) approximate the scattered field in the low frequency range as

$\begin{matrix} {{\kappa_{b}(\omega)} = {{{- \frac{1}{P_{i}}}\frac{\delta V_{cb}}{V_{b}}} = {H_{v}(\omega)}}} & ({A15}) \end{matrix}$ ${p_{s}\left( {\underline{r},\omega} \right)} = {{- {H_{p}(\omega)}}{H_{b}\left( \underline{K} \right)}\frac{e^{- {ikr}}}{r}P_{i}}$ ${H_{p}(\omega)} = \frac{K_{pL}\omega^{2}}{\omega^{2} - \omega_{0}^{2} - {i2\zeta\omega_{0}\omega}}$ $K_{pL} = {{\frac{1}{2}\rho_{t}a_{v}^{2}{LK}_{V}} = {{\frac{1}{2}{LK}_{F}} = {\frac{1}{2\omega_{s}^{2}}{LK}_{d}}}}$

For detection of the low frequency LF vibration amplitude we transmit a high frequency (HF) pulse, HF ˜10*LF or more, i.e. ω>>ω₀, as described in relation to Eqs. (12-19) and FIG. 4 . Eq. (A12) includes these HF variations of the bubble surface displacement amplitude as illustrated in 212. We can then approximate the bubble volume compressibility from Eq. (A12) for ω>>ω₀

$\begin{matrix} {{\kappa_{b}(\omega)} = {{{{H_{V}(\omega)}\frac{1}{\mu_{tr}}} \approx {{H_{d}(\omega)}\mu_{tr}} \approx {{- \frac{2{K_{d}\left( {L/a_{v}} \right)}\omega_{s}^{2}}{\omega^{2}}}\mu_{tr}}} = {{- \frac{2}{\rho a_{v}^{2}\omega^{2}}}{K_{d}\left( {L/a_{v}} \right)}}}} & ({A16}) \end{matrix}$ ${K_{d}\left( {L/a_{v}} \right)} = {\frac{1}{\pi}{\int{{dq}\frac{{qK}_{1}(q)}{K_{0}(q)}\frac{\sin^{2}{{qL}/a_{v}}}{q^{2}{L/a_{v}}}}}}$

Inserted into Eq. (13) we then get the approximation to the scattered field for the Ω_(H) frequency range as

$\begin{matrix} {{p_{s}\left( {\underline{r},\omega} \right)} \approx {{- a_{v}}{K_{pH}\left( {L/a_{v}} \right)}{H_{b}\left( \underline{K} \right)}\frac{e^{- {ikr}}}{r}P_{i}}} & ({A17}) \end{matrix}$ ${K_{pH}\left( {L/a_{v}} \right)} = {{\frac{L}{a_{v}}{K_{d}\left( {L/a_{v}} \right)}} = {\frac{1}{\pi}{\int{{dq}\frac{{qK}_{1}(q)}{K_{0}(q)}\frac{\sin^{2}{{qL}/a_{v}}}{q^{2}}}}}}$ The received HF amplitude is then proportional to Y _(r)(α_(V) ,L,K )˜α_(V) K _(pH)(L/α _(V))H _(b)( K )   (A18) With the 2^(nd) order approximation of H_(V) in Eq. (A14) the angular frequency ω_(m) for max |H_(V)(ω)| is given by

$\begin{matrix} {\zeta = {- \frac{1}{\omega_{0}{\theta^{\prime}\left( \omega_{0} \right)}}}} & ({A19}) \end{matrix}$ $\omega_{m} = {\omega_{0}\sqrt{1 - {2\zeta^{2}}}}$ $\zeta = {{\sqrt{\frac{\omega_{0}^{2} - \omega_{m}^{2}}{2\omega_{0}^{2}}}{❘{H_{V}\left( \omega_{m} \right)}❘}} = {{\frac{K_{V}}{\sqrt{\omega_{0}^{4} - \omega_{m}^{4}}}K_{V}} = {\frac{K_{d}}{\mu_{tr}} = \frac{K_{F}}{\rho a_{v}^{2}}}}}$

REFERENCES

1. Lars Hoff: “Acoustic Characterization of Contrast Agents for Medical Ultrasound Imaging” Springer Science+Business Media. 2001 

What is claimed is:
 1. A method for estimating parameters related to the surface or volume displacement vibration amplitude of at least one micro-bubble brought to vibration by an ultrasound wave of amplitude P_(L) and angular frequency ω_(L) transmitted along a LF beam with a given propagation direction along the unit vector e _(L), comprising a) transmitting towards said at least one micro-bubble at least two high frequency (HF) pulses with ω_(H)>5ω_(L) along one HF beam along the unit vector e _(H) at different times, b) directing at least one HF receive beam towards the region of the at least one micro-bubble with a direction e _(HR), the LF transmit and the HF transmit and receive beams are arranged so that they overlap at least in an observation region that includes said at least one micro-bubble, c) arranging the transmit timing of each of the transmitted HF pulses so that the pulses hits said at least one micro-bubble at different phases the LF oscillation, d) receiving at least two HF receive signals scattered from each said at least two transmitted HF pulses, and e) processing the at least two HF receive signals, to for each said at least one micro-bubble to produce an estimate of a bubble vibration parameter that gives one or both of the relative vibration of one or both of i) the bubble volume vibration relative to the resting bubble volume, and ii) a parameter related to relative bubble volume vibration wherein the HF transmit beam is directed in close to the same direction as the LF transmit beam and the LF transmit and HF transmit and HF receive beams are positioned so that the three beams overlap in a region around said at least one micro-bubble, to produce estimates of said bubble vibration parameter for each said at least one micro-bubble; a set S_(N) of N≥1 groups G_(n) of combined LF and HF pulse complexes are transmitted, each group G_(n) comprising at least two pulse complexes, where for each group the LF pulse amplitude can be zero for one pulse complex, and the LF pulse amplitude is non-zero for at least one pulse complex, and at least one of the polarity, amplitude, or phase of the LF pulse varies between pulse complexes, for each group G_(n) the HF pulse is transmitted with a timing relative t_(n) to the LF pulse so that the HF pulses hit said at least one micro-bubble with a phase relation of ω_(L)t_(n) relative to the transmitted LF oscillation, where the time t_(n) is different for different groups G_(n), the processing includes combining the HF receive signals from the scattered HF pulses from said at least two HF transmit pulses for each group G_(n) to produce a detection parameter D_(n) for group G_(n) and each said at least one micro-bubble, and for each said at least one micro-bubble, the detection parameter Dm for the group Gm that gives maximal detection parameter for all G_(n)·S_(N) is used as an estimate to obtain an estimate of the amplitude of said bubble vibration parameter.
 2. The method according to claim 1, where the HF receive signals are compensated by one or both of i) a nonlinear propagation delay, and ii) a nonlinear pulse form distortion in the processing chain to produce estimates of said bubble vibration parameter for each said at least one micro-bubble.
 3. The method according to claim 1, for an adequately dense group of micro bubbles, where a) the HF transmit beam defined by e _(H), is directed at an angle to the direction of the LF transmit beam defined by e _(L), and b) at least one HF receive beam are positioned together with the LF and HF transmit beams so that the three beams overlap at the region of said dense group of micro-bubbles, and c) transmitting a single group G_(n) comprising at least two HF pulses at a reference time t_(n), and one of ci) the LF pulse polarity is switched between at least two of the HF pulses within the group, and cii) internal timing of the pulses within the group is selected, so that at least two HF pulse complexes hits at least one bubble within said group of micro-bubbles at opposite polarities of the LF oscillation, and d) processing that includes combining the received HF signal from said at least two HF pulses to produce a detection parameter D_(n) for each bubble n in said group of micro-bubbles, and e) selecting estimate of a bubble vibration parameter the detection parameter D_(m) for the bubble n=m that has the highest value.
 4. The method according to claim 1, where said parameter related to relative bubble volume vibration is the average surface normal vibration amplitude relative to a radius of the bubble.
 5. The method according to claim 4, where the value of the average surface normal is used to set the transmit amplitude of the LF wave.
 6. The method according to claim 1, where the LF beam is focused onto in the therapeutic region.
 7. The method according to claim 1, where the HF beam is broad with width less than the LF beam in the observation region, and several multi-focused parallel HF receive beams are generated to provide spatial resolution in the estimates of said bubble vibration parameter for a group of micro-bubbles.
 8. The method according to claim 1, where scanning in 2 or 3 dimensions of one or more of i) the HF receive beam direction, and ii) the HF transmit beam direction, and iii) the LF beam direction to produce images of vibrating micro bubbles and the vibration amplitudes of the bubbles for a given transmitted amplitude of the LF wave.
 9. The method according to claim 1, where the phase of the 2^(nd) order transfer function between incident LF pressure oscillation and the amplitude bubble volume vibration is estimated from the timing t_(m) that gives the maximal detection parameter D_(m).
 10. The method according to claim 9, where the maximal detection parameter D_(m) is obtained for at least two different frequencies of the LF beam, and estimates of the following parameters of the 2^(nd) order approximation to the transfer function from incident pressure amplitude to bubble volume vibration i) the angular resonance frequency {circumflex over (ω)}₀ of the bubble ii) the absorption coefficient {circumflex over (ζ)} of the 2^(nd) order approximation to the incident LF pressure to bubble volume transfer function, iii) the gain {circumflex over (K)}_(V) of the 2^(nd) order approximation to the incident LF pressure to bubble volume transfer function.
 11. The method according to claim 10, where for a bubble filling parts of a vessel estimates of at least one of the following vibration parameters for the surrounding tissue is obtained i) the ratio of the tissue shear wave velocity c_(s) to the vessel radius a_(V) is obtained from the estimate of the angular resonance frequency {circumflex over (ω)}₀ of the bubble, ii) for a known radius of the vessel the real part of tissue shear stiffness is obtained from the estimate of the angular resonance frequency {circumflex over (ω)}₀ of the bubble, iii) for a known radius of the vessel the imaginary part of tissue shear stiffness is obtained from the estimate of the angular resonance frequency {circumflex over (ω)}₀ of the bubble.
 12. An instrument for obtaining estimates of parameters related to the surface or volume displacement vibration amplitude of at least one micro-bubble brought to vibration by an ultrasound wave of amplitude P_(L) and angular frequency ω_(L), comprising a) means for transmitting towards said at least one micro-bubble ai) at least one LF pulse along a LF transmit beam with a transmit beam direction along a unit vector e _(L), and aii) at least two high frequency (HF) pulses with ω_(H)>5ω_(L) along a HF beam along the unit vector e _(H) at different times, and iii) means for arranging the transmit timing of each of the transmitted HF pulses so that the pulses hit said at least one micro-bubble at different phases the LF oscillation, said means for transmitting include means for transmitting said HF transmit beam in close to the same direction as the LF transmit beam and further includes i) means for transmitting a set S_(N) of N≥1 groups G_(n) of combined LF and HF pulse complexes, each group G_(n) comprising at least two pulse complexes, where for each group the LF pulse amplitude can be zero for one pulse complex, and the LF pulse amplitude is non-zero for at least one pulse complex, and ii) means for varying at least one of the polarity, amplitude, or phase of the LF pulse between pulse complexes, and iii) means that for group G_(n) for sets the transmit time t_(n) of the HF pulse relative to the LF pulse so that the HF pulses hit said at least one micro-bubble with a phase relation of ω_(L)t_(n) relative to the LF pulse oscillation when the HF pulse hits the at least one micro-bubble, where the time t_(n) is different for different groups G_(n) and b) means for receiving that includes include means for directing the HF receive beam so that the LF transmit and HF transmit and HF receive beams overlap at the position of said at least one micro-bubble and further including i) means for directing at least one HF receive beam with direction e _(HR) to produce an overlap region (OR) with the LF and HF transmit beams that includes the at least one micro-bubble, and ii) means for receiving at least two HF receive signals scattered from each said at least one micro-bubble from said at least two transmitted HF pulses, and c) means for processing the at least two HF receive signals, to for each said at least one micro-bubble to produce an estimate of a bubble vibration parameter that gives one or both of the relative vibration of one or both of i) the bubble volume vibration relative to the resting bubble volume, and ii) a parameter related to relative bubble volume vibration, wherein said means for processing includes i) means for combining the scattered HF pulses from said at least two HF pulses for each group G_(n) to produce a detection parameter D_(n) for group G_(n) and each said at least one micro-bubble, and ii) means for selecting for each said at least one micro-bubble, the detection parameter D_(m) for the group G_(m), that gives maximal detection parameter for all G_(n)·S_(N), said D_(m) is used as a basis to obtain an estimate of the amplitude of said bubble vibration parameter.
 13. The instrument according to claim 12, where said processing means include means for compensating the HF receive signals by one or both of i) a nonlinear propagation delay, and ii) a nonlinear pulse form distortion in the processing chain to produce estimates of said bubble vibration parameter for each said at least one micro-bubble.
 14. The instrument according to claim 12, for an adequately dense group of micro bubbles, where a) said means for transmitting includes ai) means for directing the HF transmit beam defined by e _(H), at an angle to the direction of the LF transmit beam defined by e _(L), and aii) means for transmitting a single group G_(n) comprising at least two HF pulses at a reference time t_(n), and one of a) the LF pulse polarity is switched between at least two of the HF pulses within the group, and ii) internal timing of the pulses within the group is selected, so that at least two HF pulse complexes hits at least one bubble within said group of micro-bubbles at opposite polarities of the LF oscillation, and b) said means for receiving includes means for directing at least one HF receive beam together with the LF and HF transmit beams so that the three beams overlap at a region of said dense group of micro-bubbles, and c) said means for processing includes ci) means combining the received HF signal from said at least two HF pulses to produce a detection parameter D_(n) for each bubble n in said group of micro-bubbles, and cii) means for selecting the detection parameter D_(m) for the bubble n=m that has the highest value, and use D_(m) in the processing chain to produce and estimate of a bubble vibration parameter.
 15. The instrument according to claim 12, where the LF beam is focused onto in the therapeutic region.
 16. The instrument according to claim 12, where the HF transmit beam is broad with width less than the LF beam in the observation region, and several multi-focused parallel HF receive beams are generated to provide spatial resolution in the estimates of said bubble vibration parameter for a group of micro-bubbles.
 17. The instrument according to claim 12, where said parameter related to relative bubble volume vibration is the average surface normal vibration amplitude relative to a radius of the bubble.
 18. The instrument according to claim 17, comprising means to use the value of the average surface normal to set the transmit amplitude of the LF wave.
 19. The instrument according to claim 12, where said transmit and receive means comprises means for scanning in 2 or 3 dimensions of one or more of i) the HF receive beam direction, and ii) the HF transmit beam direction, and iii) the LF beam direction to produce images of vibrating micro bubbles and the vibration amplitudes of the bubbles for a given transmitted amplitude of the LF wave.
 20. The instrument according to claim 12, where said processing means include means to estimate the phase of the 2^(nd) order transfer function between incident LF pressure oscillation and the amplitude bubble volume vibration from the timing t_(m) that gives the maximal detection parameter D_(m).
 21. The instrument according to claim 20, where the processing means includes a) means for detection of the maximal detection parameter D_(m) for at least two different frequencies of the LF beam, and b) means for estimating the following parameters of the 2^(nd) order approximation to the transfer function from incident pressure amplitude to bubble volume vibration i) the angular resonance frequency {circumflex over (ω)}₀ for the bubble ii) the absorption coefficient {circumflex over (ζ)} of the 2^(nd) order approximation to the incident LF pressure to bubble volume transfer function, iii) the gain {circumflex over (K)}_(V) of the 2^(nd) order approximation to the incident LF pressure to bubble volume transfer function.
 22. The instrument according to claim 21, where the processing means comprises means to estimate at least one of the following vibration parameters for the tissue surrounding a bubble that is filling a part of a vessel i) the ratio of the tissue shear wave velocity c_(s) to the vessel radius a_(V) from estimates of the angular resonance frequency {circumflex over (ω)}₀ of the bubble, and ii) the real part of tissue shear stiffness for a known radius of the vessel from estimates of the angular resonance frequency {circumflex over (ω)}₀ of the bubble, and iii) the imaginary part of tissue shear stiffness for a known radius of the vessel from the estimate of the angular resonance frequency {circumflex over (ω)}₀ of the bubble.
 23. A method for estimating parameters related to the surface or volume displacement vibration amplitude of at least one micro-bubble brought to vibration by an ultrasound wave of amplitude P_(L) and angular frequency ω_(L) transmitted along a LF beam with a given propagation direction along the unit vector e _(L), comprising a) transmitting towards said at least one micro-bubble at least two high frequency (HF) pulses with ω_(H)>5ω_(L) along one HF beam along the unit vector e _(H) at different times, b) directing at least one HF receive beam towards the region of the at least one micro-bubble with a direction e _(HR), the LF transmit and the HF transmit and receive beams are arranged so that they overlap at least in an observation region that includes said at least one micro-bubble, c) arranging the transmit timing of each of the transmitted HF pulses so that the pulses hits said at least one micro-bubble at different phases the LF oscillation, d) receiving at least two HF receive signals scattered from each said at least two transmitted HF pulses, and e) processing the at least two HF receive signals, to for each said at least one micro-bubble to produce an estimate of a bubble vibration parameter that gives one or both of the relative vibration of one or both of i) the bubble volume vibration relative to the resting bubble volume, and ii) a parameter related to relative bubble volume vibration wherein the HF transmit beam is directed in close to the same direction as the LF transmit beam and the LF transmit and HF transmit and HF receive beams are positioned so that the three beams overlap in a region around said at least one micro-bubble, to produce estimates of said bubble vibration parameter for each said at least one micro-bubble; wherein a long LF pulse of duration T_(L) is transmitted towards said at least one micro-bubble, where T_(L) is longer than T_(R), preferably longer than N*TR, where T_(R) is the propagation time from HF transmit transducer to said at least one micro-bubble, and a set S_(N) of N≥1 groups G_(n), comprising two HF pulses are transmitted at a reference time t_(n) for each group, and internal timing within each group is selected so that at least two pulse complexes hits the bubble at opposite polarities of the LF oscillation, and the processing includes combining the scattered HF pulses from said at least two HF pulses for each group G_(n) to produce a detection parameter D_(n) for group G_(n) and each said at least one micro-bubble, and for each said at least one micro-bubble, the detection parameter Dm for the group G_(m), that gives maximal detection parameter for all G_(n)·S_(N) is used as an estimate to obtain an estimate of the amplitude of said bubble vibration parameter.
 24. The method according to claim 23, where the HF receive signals are compensated by one or both of i) a nonlinear propagation delay, and ii) a nonlinear pulse form distortion in the processing chain to produce estimates of said bubble vibration parameter for each said at least one micro-bubble.
 25. The method according to claim 23, where said parameter related to relative bubble volume vibration is the average surface normal vibration amplitude relative to a radius of the bubble.
 26. The method according to claim 25, where the value of the average surface normal is used to set the transmit amplitude of the LF wave.
 27. The method according to claim 23, where the LF beam is focused onto in the therapeutic region.
 28. The method according to claim 23, where the HF beam is broad with width less than the LF beam in the observation region, and several multi-focused parallel HF receive beams are generated to provide spatial resolution in the estimates of said bubble vibration parameter for a group of micro-bubbles.
 29. The method according to claim 23, where scanning in 2 or 3 dimensions of one or more of i) the HF receive beam direction, and ii) the HF transmit beam direction, and iii) the LF beam direction to produce images of vibrating micro bubbles and the vibration amplitudes of the bubbles for a given transmitted amplitude of the LF wave.
 30. An instrument for obtaining estimates of parameters related to the surface or volume displacement vibration amplitude of at least one micro-bubble brought to vibration by an ultrasound wave of amplitude PL and angular frequency ω_(L), comprising a) means for transmitting towards said at least one micro-bubble ai) at least one LF pulse along a LF transmit beam with a transmit beam direction along a unit vector e _(L), said LF pulse comprising a long LF pulse of duration T_(L) towards said at least one micro-bubble, where T_(L) is long, preferably longer than NT_(R) where T_(R) is the propagation time from HF transmit transducer to said at least one micro-bubble, and a set S_(N) of N≥1 groups G_(n) comprising two HF pulses transmitted at a reference time t_(n) for each group, and selecting internal timing within each group so that at least two pulse complexes hits the bubble at opposite polarities of the LF oscillation and aii) at least two high frequency (HF) pulses with ω_(H)>5 ω_(L) along a HF beam along the unit vector e _(H) at different times, and aiii) means for arranging the transmit timing of each of the transmitted HF pulses so that the pulses hit said at least one micro-bubble at different phases the LF oscillation, wherein said means for transmitting include means for transmitting said HF transmit beam in close to the same direction as the LF transmit beam, and b) means for receiving that includes i) means for directing at least one HF receive beam with direction e _(HR) to produce an overlap region (OR) with the LF and HF transmit beams that includes the at least one micro-bubble, and ii) means for receiving at least two HF receive signals scattered from each said at least one micro-bubble from said at least two transmitted HF pulses, and wherein said means for receiving include means for directing the HF receive beam so that the LF transmit and HF transmit and HF receive beams overlap at the position of said at least one micro-bubble; and c) means for processing the at least two HF receive signals, to for each said at least one micro-bubble to produce an estimate of a bubble vibration parameter that gives one or both of the relative vibration of one or both of i) the bubble volume vibration relative to the resting bubble volume, and ii) a parameter related to relative bubble volume vibration, wherein said means for processing includes means for combining the received scattered HF pulses from said at least two HF transmit pulses for each group Gh_(n) to produce a detection parameter D_(n) for group G_(n) and each said at least one micro-bubble, and means for selecting for each said at least one micro-bubble the detection parameter D_(m) for group G_(m), that gives maximal detection parameter for all G_(n)·S_(N), Said G_(m) is used to obtain an estimate of the amplitude of said bubble vibration parameter for each said at least one micro-bubble.
 31. The instrument according to claim 30, where said processing means include means for compensating the HF receive signals by one or both of i) a nonlinear propagation delay, and ii) a nonlinear pulse form distortion in the processing chain to produce estimates of said bubble vibration parameter for each said at least one micro-bubble.
 32. The instrument according to claim 15, where the LF beam is focused onto in the therapeutic region.
 33. The instrument according to claim 16, where the HF transmit beam is broad with width less than the LF beam in the observation region, and several multi-focused parallel HF receive beams are generated to provide spatial resolution in the estimates of said bubble vibration parameter for a group of micro-bubbles.
 34. The instrument according to claim 17, where said parameter related to relative bubble volume vibration is the average surface normal vibration amplitude relative to a radius of the bubble.
 35. The instrument according to claim 34, comprising means to use the value of the average surface normal to set the transmit amplitude of the LF wave.
 36. The instrument according to claim 30, where said transmit and receive means comprises means for scanning in 2 or 3 dimensions of one or more of i) the HF receive beam direction, and ii) the HF transmit beam direction, and iii) the LF beam direction to produce images of vibrating micro bubbles and the vibration amplitudes of the bubbles for a given transmitted amplitude of the LF wave. 